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We construct improved quantum Monte Carlo estimators for the spherically- and system-averaged 
electron pair density (i.e. the probability density of finding two electrons separated by a relative 
distance u), also known as the spherically-averaged electron position intracule density I{u), using 
the general zero-variance zero-bias principle for observables, introduced by Assaraf and Caffarel. 
The calculation of I{u) is made vastly more efficient by replacing the average of the local delta- 
function operator by the average of a smooth non-local operator that has several orders of magnitude 
smaller variance. These new estimators also reduce the systematic error (or bias) of the intracule 
density due to the approximate trial wave function. Used in combination with the optimization of 
an increasing number of parameters in trial Jastrow-Slater wave functions, they allow one to obtain 
well converged correlated intracule densities for atoms and molecules. These ideas can be applied 
to calculating any pair-correlation function in classical or quantum Monte Carlo calculations. 



I. INTRODUCTION 

Two-electron distribution functions occupy an impor- 
tant place in electronic structure theory between the sim- 
plicity of one-electron densities and the complexity of the 
many-electron wave function. In particular, the system- 
averaged electron pair density, the probability density of 
finding two electrons separated by the relative position 
vector u, also known as the electron position intracule 
density /(u), plays an important role in qualitative and 
quantitative descriptions of electronic systems. Position 
intracule densities have been extensively used to analyze 
shell structure, electron correlation, Hund's rules and 
chemical bonding (see, e.g., Refs. 1-31). Density func- 
tional theory-like approaches have been proposed based 
on the position intracule density [32-36] or on the closely- 
related Wigner intracule density [37-39] . 

Position intracule densities have been extracted from 
experimental X-ray scattering intensities for small atoms 
and molecules [40-42]. They have been calculated in 
the Hartree-Fock (HF) approximation for systems rang- 
ing from small atoms to large molecules [7, 12, 16, 
19, 21, 27, 43-46]. Calculations using common quan- 
tum chemistry correlated methods such as second-order 
M0ller-Plesset perturbation theory, multi-configurational 
self-consistent-field (MCSCF) and configuration interac- 



tion approaches have been limited to atoms and small 
molecules [6, 11, 14, 18, 24, 28, 47, 48]. Very accurate cal- 
culations using Hylleraas-type explicitly-correlated wave 
functions have been done only for the helium and lithium 
isoelectronic series [1, 5, 8, 15, 17, 49-53]. Variational 
Monte Carlo (VMC) calculations using Jastrow-Slater 
wave functions have been used to compute correlated po- 
sition intracule densities for atoms from helium to neon 
and some of their isoelectronic series [22, 23, 25, 26, 29- 
31, 52, 54, 55]. In this paper, we show that the calcula- 
tion of position intracule densities using quantum Monte 
Carlo (QMC) methods can be made much more accurate 
and efficient, opening new possibilities of investigation. 

The position intracule density associated with an 
TV-electron (real) wave function ^'(R), where R — 
(ri,r2, ...jTn) is the 3-/V— dimensional vector of electron 
coordinates (ignoring spin for now), is defined as the 
quantum-mechanical average of the delta-function oper- 
ator (5(r,y — u) 



(1) 



where = Vj — and the sum is over all electron pairs, 
and its spherical average is 
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Among the most important properties of I{u) are the 
normalization sum rule (giving the total number of elec- 
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tron pairs) 
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the electron-electron cusp condition [56] 
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and, for finite systems, the exponential decay at large 
u determined by the spherically averaged one-electron 
density n{u) evaluated at the distance u from the chosen 
origin 
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where / is the vertical ionization energy (see Refs. 57- 
59). The moments of the radial intracule density 

Mfe = Jp°° d.uA'mP'^^ I{u) are related to physical observ- 
ables [60], in particular M_i is just the electron-electron 
Coulomb interaction energy 
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In standard correlated methods based on an expansion 
of the wave function in Slater determinants, the impor- 
tant short-range part of the position intracule density 
converges very slowly with the one-electron and many- 
electron basis. The advantage of employing QMC meth- 
ods [61] lies in the possibility of using compact, explicitly- 
correlated wave functions which arc able to describe 
properly the short-range part of I{u). However, the prob- 
lem in this approach is that one calculates the average of 
a delta-function operator which has an infinite variance. 
As for other probability densities, the standard procedure 
in QMC approaches simply consists of counting the num- 
ber of electron pairs separated by the distance u within 
Au encountered in the Monte Carlo run. More precisely, 
e.g. in VMC calculations, I{u) is estimated as the sta- 
tistical average 



(7) 



over the trial wave function density ^'(R)^, of the local 
"histogram" estimator /]J"'*°(m,R) 
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where l[a, &[(''') = 1 if r G [a,b[ and otherwise. Here 

and in the following, (/(R))^2 = (1/M) E^i /(Rfe) 
designates the statistical average of the local quantity 
/(R) over M configurations Rk sampled from 5'(R)^. 
For M ^ 0, this histogram estimator has a finite but 
large variance for small u or small Au. Although the 
use of importance sampling can help decrease the vari- 
ance [22, 62, 63], the calculation of position intracule 



densities in QMC remains very inefficient: very long runs 
have to be performed to reach an acceptably small sta- 
tistical uncertainty. Moreover, the histogram estimator 
has a discretization error: even in the limit of an infinite 
sample M ^ oo, J'"'**°(u) remains only an approximation 
of first order in Au to the position intracule density I{u) 
of the wave function \1/ (R) . Note however that it is pos- 
sible to greatly reduce the discretization error by choos- 
ing a flexible analytic form for I{u) that obeys known 
conditions (such as the cusp condition of Eq. (4)) and 
fitting in each interval Aw the integral of I{u) (instead 
of I{u)) to the computed data. This approach has been 
used to calculate radial electron densities for atoms [64] 
and spherically-averaged intracule densities [52] but the 
method is not as effective for multidimensional densities. 

In this work, we develop improved QMC estimators 
for the spherically averaged position intracule density 
based on the zero-variance zero-bias (ZVZB) principle 
for observables introduced by Assaraf and Caffarel for 
calculations of electronic forces [65-67], which has been 
recently applied to calculations of one-electron densi- 
ties [68]. These new ZVZB estimators of /(u) have vari- 
ances several orders of magnitude smaller than the ones 
obtained with the standard estimator of Eq. (8) and thus 
dramatically increase the efficiency of these calculations. 
Moreover, these estimators do not suffer from any dis- 
cretization error and can even reduce the systematic error 
due to the approximate trial wave function. Like related 
techniques proposed for calculations of one-electron or 
two-electron densities using deterministic ab initio meth- 
ods [52, 69-83], these estimators replace the average of 
the local delta-function operator by the average of an 
operator which is nonlocal in real space. They can be 
viewed as a generalization of the improved QMC estima- 
tors previously proposed for computations of averages of 
probability densities at particle coalescences [62, 84 86]. 

Moreover, we make use in this work of the recently- 
developed linear energy minimization method [87, 88] to 
finely optimize the Jastrow parameters, the configuration 
state functions (CSFs) coefficients and the orbital coef- 
ficients of our trial Jastrow-Slater wave functions. Op- 
timization of the determinantal part of the wave func- 
tion with an increasing number of CSFs allows us to ob- 
tain a systematic improvement of wave functions. This 
provides a practical route for calculating intracule densi- 
ties in VMC and fixed-node (FN) diffusion Monte Carlo 
(DMC) with progressively smaller systematic errors. 

The paper is organized as follows. In Sec. II, we review 
the principle of zero- variance zero-bias improved estima- 
tors for an arbitrary observable in QMC. In Sec. Ill, we 
give improved estimators for the case of the position in- 
tracule density. Sec. IV contains computational details 
of the calculations, and Sec. V discusses results for the 
He and C atoms and for the C2 and N2 molecules to il- 
lustrate the technique. Sec. VI contains our conclusions. 

Hartree atomic units are used throughout this work. 
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II. ZERO- VARIANCE ZERO-BIAS IMPROVED 
ESTIMATORS 

Wc review here the principle of the zero- variance zero- 
bias improved QMC estimators for an arbitrary observ- 
able that does not commute with the Hamiltonian, de- 
veloped in Refs. 65-67. Throughout this section, all the 
averages (• • • are considered in the limit of an infinite 
Monte Carlo (MC) sample M ^ oo. 



A. Estimators in variational Monte Carlo 

In VMC, the exact energy Eq = (*o|-ff |*o)/(*o|^'o) 
of some exact eigcnfunction 4*0 of the electronic Hamil- 
tonian H is estimated by the statistical average of the 
local energy El{R) = (R|^|*)/*(R) 



(9) 



using an approximate trial wave function 'I'(R). The 
systematic error (or bias) of this estimator SE = E — 

Eq and its variance cr^ (El) = ((-El(R) — -E')^)*2, whose 
square root is proportional to the statistical uncertainty, 
both vanish quadratically as a function of the error in 
the trial wave function = l^t — "Jfol (where |---| 

designates the Hilbert space norm) 



SE = 0{\5'i/\^), 



{El) = Om\% 



(10a) 
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which is referred to as the quadratic zero- variance zero- 
bias property of the local energy. This is easily shown by 
writing "if = "ifo + S'i^ in the expressions of the average 
and the variance, and expanding to second-order in S^. 

Similarly, the exact expectation value of an arbitrary 
observable Oo = (*o|0|*o)/(*o|*o) can be estimated by 
the statistical average of the local observable Oi,(R) = 
(R|d|*)/*(R) 



O = (Ol(R))^2 



(11) 



but, as ^0 is generally not an eigenstate of O, the sys- 
tematic error of this estimator 50 = O — Oo vanishes 
only linearly with |(5^|, while its variance cr^ {Ol) = 
((Ol(R) — O) generally does not even vanish with 
\S^\ 



60 = 0{\6^\), 



(Ol) = 0{1), 



(12a) 



(12b) 



which often makes the calculations of observables inac- 
curate and inefficient. 

However, Assaraf and Caffarel [67] have pointed out 
that the quadratic zero- variance zero-bias property of the 
energy can be extended to an arbitrary observable by 



expressing it as an energy derivative. This is based on 
the Hellmann-Feynman theorem which states that, if one 
defines the A-dependent Hamiltonian 



H^ = H + X0, 



(13) 



with some associated exact A-dependent eigcnfunction 
*^ = *o + A*'o + --- , (14) 

where ^'q = (d^'g /dA) and corresponding exact A- 

dependent energy E^ = then the 

exact value of the observable is given by the deriva- 
tive of the energy with respect to A at A = 0: Oo = 
[cIEq / dX) Introducing an approximate A-dependent 
trial wave function 



= * + A*' 



(15) 



where ^' = [d^^/dX) the A-dependent energy can 
be estimated as the statistical average of the A-dependent 
local energy E^(R) = (R|fl"^|*^)/*^(R) 



(16) 



over the probability density \E''*'(R)^, and the Hellmann- 
Feynman theorem suggests now to define a zero-variance 
zero-bias (ZVZB) estimator for the observable as the 
derivative of E^ with respect to A at A = 

qZvzb ^ (or^^iK))^, 



dE^ 
~dX 



A=0 



= (04R))^2 + (A0r(R))* 



+ (AO|^R))^,, 
where the zero-variance (ZV) contribution 
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(RI.H'I*') 
*'(R) 



El{-R) 



¥(Ry' 



(17) 



(18) 



does not contribute to the average, (AO^^(R))^2 = 
(in the limit of an infinite MC sample M oo, from 
the Hermiticity of the Hamiltonian) but can decrease the 
variance, and the zero-bias (ZB) contribution 



AOf (R)=2[i?i(R)-i;] 



^^(R) 
*(R) ' 



(19) 



usually has little effect on the variance but can decrease 
the systematic error. The average of this ZB term van- 
ishes if the trial wave function ^ is an exact eigenstate or, 
more generally, if the energy E is stationary with respect 
to an infinitesimal variation of the trial wave function 
along the derivative ^'i ^ ^ ^ -\- e^' . For the special 
case of the calculation of the derivatives of the energy 
with respect to the nuclear coordinates, the average of 
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the ZB term is known in the electronic structure litera- 
ture as the Pulay contribution [89]. 

If one wants simply to compute an observable for a 
given trial wave function without changing the average, 
the ZB correction has to be omitted. In this case, the 
ZV estimator Ol^{R) = Ol(R) + AO|^(R) has a sys- 
tematic error that still vanishes linearly with jjvl'l but a 
variance that now vanishes quadratically with and 
the error in its derivative |i5^'| = |^'' — 

50zv ^ 0{\d^\), (20a) 

(Ofv^ ^ 0(|5'J'|2 + |5^-'|2 + i^^-i i^^'i). (20b) 

If one also includes the ZB correction, then both the sys- 
tematic error and the variance of the ZVZB estimator 
0|^^^(R) vanish quadratically with {S'l'l and 

<50^^2^ = 0(1(5*12 + |(54'||,5*'|), (21a) 

^2 (oZVZB^ ^ OdJ^-p + |55.'|2 + l^^l |<5^'|). (21b) 

This is easily shown by writing ^I* = 5'o + f^^* and 
^' = ^'q + S^' in the expressions of the average and the 
variance, and expanding with respect to and 6"^' . One 
can verify that the exact cancellation of the dominant 
contributions to the systematic error and the variance of 
Eqs. (12a) and (12b) by the ZV and ZB corrections is 
a direct consequence of the relationship between the ex- 
act wave function and its derivative in first-order 
perturbation theory: {H - Eo)\^'q) = -(6 - Oo)|^'o)• 
Thus, the ZVZB estimators permit in principle accu- 
rate and efficient calculations for observables, provided 
that a good approximation for the derivative \E'q is avail- 
able. We note that optimizing the parameters p in the 
trial wave function 'J'(p) by energy minimization presents 
the advantage of decreasing the sensitivity of the system- 
atic error to the quality of the approximate derivative ^' , 
since in this case only needs to be a good approxima- 
tion of in the orthogonal complement of the space 
spanned by the wave function derivatives with respect 
to the parameters at the optimal parameter values Poptj 
i.e. = (94'(p)/9pi)p^p^^^. Indeed, the contribution to 
the average of the ZB term of any component of ^' along 
^'i is proportional to the energy gradient with respect to 
parameter pi and thus vanishes at the stationary point: 

(2[i?i(R)-i;]*i(R)/*(R)),„ = (aii;/aK)p=p,,, = o. 

This point is at the origin of the improved accuracy in 
the calculations of forces in QMC in Refs. 90, 91. 

B. Estimators in diffusion Monte Carlo 

DMC calculations within the FN approximation go 
beyond VMC calculations by replacing the VMC dis- 
tribution ^'(R)2 by the more accurate mixed FN-DMC 
distribution \1/fn(R)*(R) in statistical averages, i.e. 



(• • •)^2 (• • O^FN*- the DMC variances have the 
same lowest order behaviors with respect to the error in 
the trial wave function as in VMC. The systematic error 
of the energy now vanishes as the product of the error 
in the trial wave function |(5^'| and the error in the FN 
wave function |^^'fn| = |*fn — ^o| 

SE = 0{\5^\ |(5*fn|). (22) 

For an arbitrary observable, the systematic error still 
vanishes only linearly with \d^\ or |^'J'fn| 

SO = 0{\6^\ + |5*fn|), (23) 

but the use of the hybrid estimator "2DMC - VMC", 
Ohybrid = 2(Oz,(R))^p^^ - {Ol{'R))^2, allows one to 
remove the dominant linear term \S^\ 

(50hybrid = 0{\d^F^\ + |,5*|2 + |(5*fn|' + |(5*fn|). 

(24) 

The same improved estimators defined in the previous 
section can also be used straightforwardly in FN-DMC 
calculations. Note that, in this case, the average of the 
ZV term of Eq. (18) no longer vanishes on an infinite MC 
sample. Nevertheless, this term gives the same order of 
reduction of the variance as in VMC. The systematic 
error and the variance are still given by Eqs. (20a) and 
(20b), though the prefactors can be different in VMC 
and DMC. We note that it is possible to define a new ZV 
correction of vanishing average in FN-DMC [66] but this 
does not change the leading order of either the systematic 
error or the variance. 

Adding the ZB term of Eq. (19), the systematic error 
of the ZVZB estimator vanishes quadratically 

^qZvzb ^ Oi\S^\^ + \S^\\S^'\ + \S^\\6<bFN\ 

+ \S-9'\ |<5*fn|), (25) 

with no linear term |(5\1/fn| in contrast to the hybrid 
estimator of Eq. (24). The behavior of the variance 
of the ZVZB estimator in FN-DMC is still given by 
Eq. (21b). One can also define a hybrid ZVZB estimator, 
Oh^b^^fd = 2{01^^^{K))^^^^ - (Or^^(R))^., whose 
systematic error also vanishes quadratically as 

SOl^h^a = 0(|<5*|2 + |^*||5fFN| + |^*'||^*FN|), 

(26) 

with no term in \6^\ \S^'\ in contrast to Eq. (25). 

C. Expressions of the estimators for local 
Hamiltonians 

We give now more explicit expressions of the estimators 
in terms of the convenient logarithmic derivative Q(R) = 
*'(R)/*(R). 

In the case of local Hamiltonians (R|_ff|R') = 
H{R)6(R-R') where iJ(R) = -(1/2) Ef=i V2,+y(R) 
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contains kinetic and potential energy contributions, the 
ZV term takes the form 



1 ^ 

AOr(R) = -2E[^r.Q(R) + 2vfe(R).V: 



fe=l 



(27) 



where Vfe(R) = Vrfc*(R)/*(R) is the drift velocity of 
the trial wave function. The ZB term is simply 



AOfS(R) =2A£;i(R)g(R), 



(28) 



where AEl{R) = El(R) — E. For each observable, a 
specific form for Q(R) has to be chosen. In principle, 
one could optimize (3(R) for eac:li system by minimizing 
the variance of the ZVZB estimator. Besides (5(R), the 
ZV and ZB terms require only the evaluation of the drift 
velocity and the local energy of the trial wave function 
that are already available in QMC programs. 

Clearly, a nonlocal operator y(R, R') in the Hamil- 
tonian, e.g. a nonlocal pseudopotential, would yield an 
additional term in the simplified expression of the ZV 
correction of Eq. (27) . Although the inclusion of this ad- 
ditional term is required for the variance to rigorously 
vanish quadratically as in Eq. (20b), in practice the es- 
timator using the simple form of Eq. (27) can already 
achieve a considerable variance reduction. 



III. ZERO- VARIANCE ZERO-BIAS IMPROVED 

ESTIMATORS FOR THE SPHERICALLY 
AVERAGED POSITION INTRACULE DENSITY 

In the case of the spherically averaged position intrac- 
ule density, the local observable II {u, R) is 

hiu,n) = lY: I ^^(r..-")' (29) 

which has an infinite variance. We will give two possible 
choices for Q(R) which define good improved estimators. 



A. First improved estimator 

The minimal choice for Q(R) that cancels the infinite 
variance of II {u, R) is 



dQ,u 



477 |ry-u| 



(30) 



which gives the following ZV contribution 



Sir J Att 

-F2Vi(R)- Vr 



r,,- - u 



(31) 



Using the identity .l/|rjj — u| 



-47r5(rj 



u), it 



is seen that the first term in A/^^^(u, R) exactly can- 
cels the delta function in Il{u, R), and the ZV estimator 
Il^^{u, R) = Il{u, R) + All'^^u, R) is simply 



tZVI 



(w,R) 



47r ^ 



.(R) 



dflu 

4-77 



u 



r - ■ 



The form of (5i(u, R) of Eq. (30) could also have been 
guessed from the behavior at small rjj of the first-order 
solution \E''(R) of the perturbation theory with respect to 
(5(ry — u) (see Refs. 72, 92). In contrast to the histogram 
estimator of Eq. (8) where only electron pairs with rela- 
tive distance in the small interval [u — Au/2, u+Au/2[ 
contribute to I{u), for the ZV estimator of Eq. (32) all the 
electron pairs with relative distance Vij > u contribute 
to I{u). For u 0, the variance of this ZV estimator 
is finite and much smaller than the variance of the his- 
togram estimator as shown in Sec. V. For u = 0, this ZV 
estimator reduces to the estimator proposed in Ref. 62 
to compute one-electron densities at the nucleus and ap- 
plied in Ref. 22 for calculations of intracule densities at 
the coalescence point. 

The ZB correction corresponding to (5i(w, R) of 
Eq. (30) is 



A/|«1(m,R) = - 



AEi(R) 



E 



An 

AEl(R) ^ iMoMi^ij) 



An 



u 



47r 



■E 



1 



u,-\-oo 



(33) 



The resulting ZVZB estimator 



ZVIZBI 



(m,R) 



/f'^^(u,R) + A/|^1(m, R) allows one to reduce the sys- 
tematic error due to the trial wave function as shown 
in Sec. V. Note that, in contrast to the histogram esti- 
mator, these ZV and ZVZB improved estimators do not 
suffer from any discretization error: they are estimators 
of I{u) for an infinitely precise value of u. Both the ZV 
estimator of Eq. (32) and its ZB correction of Eq. (33) 
are very simple to program and very fast to compute on 
a grid over u. 

Like the histogram estimator, the ZV estimator of 
Eq. (32) gives simply zero for distances u beyond the 
largest electron-electron distance rij encountered in the 
Monte Carlo run. Also, because the form of Qi(m, R) 
of Eq. (30) has been chosen only to cure the deficiency 
of the histogram estimator at small u, the ZB correc- 
tion of Eq. (33) actually makes the histogram estimator 
worse at large u. The ZB correction decays only as 1 /u 
as u — > 00 instead of the correct exponential decay and 
consequently gives large variances at large u. Thus, it is 
inaccurate to extract the long-range behavior of the in- 
tracule density using the ZVZB estimator lf^^'^^^{u, R). 
The second choice for Q{R) presented next remedies this 
limitation. 
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B. Second improved estimator 



IV. COMPUTATIONAL DETAILS 



A more general choice for (5(R) that cancels the infi- 
nite variance of /^(u, R) is 



Att \rij - u| ' 



(34) 



where g{Yij,n) is some function satisfying g{vij,\i) = 1 
for Tij = u. For example, using 



g{Yij,M) = e 



^ -C|r.,-u| 



(35) 



with C > 0, has the advantage of ensuring a correct ex- 
ponential decay of the estimator at infinity for a finite 
system. 

The corresponding ZV estimator is 



Stt 

1 



1 



2Vi(R)- Vr 



-C|ry-u| 



ri,- - u 



ij 

e 



An '-^ L r' 



and its ZB correction is 

A£;i,(R) ^ f dfl^ e-^l''*^-"! 



(36) 



A/r^(«,R) 



Ait 



■E 



47r |rjj — u| 



(37) 



where 



A{rij,u) = 



e-C\rij-u\ _ g-f(r-y+M) 

2CrijU 



(38) 



We have chosen to illustrate the efficiency and ac- 
curacy of the new improved QMC estimators by cal- 
culating the intracule density of the He and C atoms 
and the C2 and N2 molecules in their ground states. 

The molecules are considered at their experimental ge- 
ometries (dc-c = 2.3481 Bohr [93] and (Jn-n = 2.075 
Bohr [94]). 

We start by generating a standard ab initio 
wave function using the quantum chemistry program 
GAMESS [95], typically a HF wave function or a MC- 
SCF wave function with a complete active space (CAS) 
generated by distributing n valence electrons in m va- 
lence orbitals [CAS(n,m)]. For each system considered, 
we choose the one-electron Slater basis so as to ensure 
that the HF intracule density is reasonably converged 
with respect to the basis. For the He atom, we use the 
basis of Clementi and Roctti [96], with exponents reop- 
timized at the HF level by Koga et al. [97]. For the C 
atom, we use the CVBl basis of Ema et al. [98]. For the 
C2 and N2 molecules, we use the CVB2 basis of the same 
authors. In GAMESS, each Slater function is actually 
approximated by a fit to 14 Gaussian functions [99-101] 

This standard ah initio wave function is then multi- 
plied by a Jastrow factor, imposing the electron-electron 
cusp condition, but with essentially all other free pa- 
rameters chosen to be zero to form our starting trial 
Jastrow-Slater wave function, and QMC calculations are 
performed with the program CHAMP [102] or QMC- 
MOL [103] using the true Slater basis set rather than 
its Gaussian expansion. The Jastrow parameters, the 
orbital coefficients and the configuration state func- 
tions (CSFs) coefficients of this trial wave function are 
optimized in VMC using the very efficient, recently- 
developed linear energy minimization method [87, 88] 
and an accelerated Metropolis algorithm [104, 105]. Once 
the trial wave function has been optimized, we com- 
pute the intracule density in VMC, and in DMC using 
the fixed-node and the short-time approximations (see, 
e.g., Refs. 106-109). We use an imaginary time-step of 
T = 0.01 hartree"^ in an efficient DMC algorithm featur- 
ing very small time-step errors [110]. 



RESULTS AND DISCUSSION 



B{nj,u) 



sgn(rij - u)e' 



-C\rij-u\ 



,(39) 



and where sgn is the sign function. These refined ZV esti- 
mator ll'^^u,'R.) and ZVZB estimator if^'^^^^(u.'R) = 
If^^ {u, R) + A/f ^2 {u, R) now both decay correctly expo- 
nentially as u — > 00. The constant C, is chosen according 
to Eq. (5), i.e. C, = 2^/21 where / is an estimate of the 
vertical ionization energy. On the other hand, the com- 
putational cost is larger per Monte Carlo step than that 
of the simpler estimators of the previous subsection. 



A. He atom 

The improvement due to the ZV and ZVZB estimators 
is illustrated for the simple case of He atom. 

Figure 1 shows the spherically averaged intracule den- 
sity I{u) as a function of the electron-electron distance 
u, calculated in VMC using the standard histogram esti- 
mator of Eq. (8), the ZVl improved estimator of Eq. (32) 
and the ZVIZBI improved estimator of Eq. (33), employ- 
ing only 100 000 MC configurations sampled from a HF 
trial wave function (without a Jastrow factor). In this 



7 



0.25 



- 0.15 



histogram estimator witii HF wave function 
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FIG. 1: Spherically averaged position intracule density I{u) 
as a function of the electron-electron distance u for the He 
atom. The intracule densities calculated in VMC with the his- 
togram, ZVl and ZVIZBI estimators using only 100 000 MC 
configurations sampled from a HF trial wave function (with- 
out a Jastrow factor) are compared. For the histogram esti- 
mator, a grid step of Au = 0.005 has been used. The intrac- 
ule density obtained with an accurate, explicitly-correlated 
Hylleraas-type wave function is also shown. 



ZVIZBI estimator witli HF wave function 
ZV2ZB2 estimator witii HF wave function 
accurate intracule 




-0.0002 



3.2 



3.4 3.6 
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FIG. 2: Long-range tail of the spherically averaged position 
intracule density I{u) of the He atom. The intracule densities 
calculated in VMC with the ZVIZBI and ZV2ZB2 (with ( = 
2V2I ~ 2.7) estimators using only 100 000 MC configurations 
sampled from a HF trial wave function (without a Jastrow 
factor) are compared. The intracule density obtained with an 
accurate, explicitly-correlated Hylleraas-like wave function is 
also shown. 



subsection, we intentionally use a poor trial wave func- 
tion to demonstrate the large reduction in the bias re- 
sulting from the ZVZB estimator. The calculations take 
a few seconds on a present-day single-processor personal 
computer. The intracule density obtained with an accu- 
rate, explicitly-correlated Hylleraas-type wave function 
with 491 terms [52, 111, 112] is also shown as a reference. 
The statistical uncertainty obtained with the histogram 
estimator is very large, especially at small u, making it 
impossible to extract the on-top intracule density /(O). 
The ZVl estimator spectacularly reduces the statistical 
uncertainty which becomes invisible at the scale of the 
plot, its variance being smaller by nearly 4 orders of mag- 
nitude for small and large u and by about 2 orders of 
magnitude for intermediate u. The ZVIZBI estimator 
in turn spectacularly reduces the systematic error due to 
the use of a HF trial wave function, the obtained intrac- 
ule density agreeing well with the accurate reference. In 
particular, the maximum at small u caused by the cusp at 
M = is accurately reproduced. On the other hand, the 
use of the ZVIZBI estimator increases slightly the vari- 
ance at small u and more importantly at large u where 
the ZBl correction of Eq. (33) decays too slowly as 1/u. 

The behavior of the estimators in the long-range tail 
of I{u) is explored in detail in Fig. 2 which compares the 
ZVIZBI and ZV2ZB2 estimators of Eqs. (36) and (37) 
for M > 3. As mentioned in Sec. HIB, the relative statis- 
tical uncertainty obtained with the ZVIZBI is very large 
(about 200% at u = 3.5), whereas the ZV2ZB2 estimator, 
which has the correct exponential decay at large u, spec- 
tacularly reduces the statistical uncertainty, its variance 
being smaller by about 3 orders of magnitude at u = 3.5. 



The improvement is even more dramatic for larger u. 

We now discuss the scaling of the intracule density 
and its statistical uncertainty with respect to the nuclear 
charge Z. For this purpose, we have calculated the intrac- 
ule densities of the elements of the He isoelectronic series 
from Z = 2 (He) to Z = 20 (Ari8+) using Jastrow-Slater 
wave functions. We found that, in the relevant range of 
u (it < 5), the intracule density roughly obeys the scal- 
ing law I{u/Z) ~ Z'"- with n « 3, in agreement with the 
value predicted by a simple hydrogenic model [113]. The 
statistical uncertainty of I{u/Z) computed with the ZVl 
estimator roughly obeys the same law, so that the rela- 
tive statistical accuracy on I{u/Z) does not deteriorate 
with Z in the series (at least up to Z = 20). 

Finally, we note in passing that Fig. 2 also reveals an- 
other advantage of our improved estimators: because the 
estimator of I{u) is strongly correlated with the estima- 
tor of I{u + Au), as obvious for instance in Eq. (32), the 
obtained intracule density curves are always very smooth, 
even at a scale much smaller than the statistical uncer- 
tainty. This makes the calculated intracule densities di- 
rectly suitable for subsequent manipulations. 

B. C atom 

We continue the discussion of the improved estimators 
for the more interesting case of the C atom. 

The repercussion of the shell structure of the system 
on the intracule is apparent in Fig. 3 which shows the in- 
tegrand of the electron-electron Coulomb interaction en- 
ergy over the electron-electron distance tt, i.e. A'kuI{u) 
[see Eq. (6)]. The ZVZB2 estimator is used with an ac- 
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FIG. 3: Integrand of the electron-electron Coulomb inter- 
action energy 47ru/(ii) over the electron-electron distance u 
for the C atom, calculated in VMC using the ZVZB2 (with 
C, = 2y/2J « 1.82 [114]) improved estimator and a fully opti- 
mized Jastrow-Slater CAS (4,4) trial wave function. The same 
spin contribution 4-!tuI^^{u) and the opposite spin contribu- 
tion 47rtt/°^(u) are also shown. 

curate Jastrow-Slater CAS (4,4) wave function in which 
the Jastrow, CSF and orbital parameters have been si- 
multaneously optimized. Note however that, at the scale 
of the plot, the HF intracule 47ru/HF('u) would look the 
same. The curve displays two maxima: the maximum 
at short distance (m ^ 0.2) corresponds essentially to 
the K-K electron pair and the maximum at longer dis- 
tance {u Ri 1.0) corresponds essentially to the K-L and 
L-L electron pairs. Fig. 3 also shows the same spin (SS) 
contribution 

J^'5(u) = /"(u)+/ii(u), (40) 

and opposite spin (OS) contribution 

/°S(u) = /Ti(u) + /iT(y)^ (41) 

where I'"^ (u) is the spherical average of the spin-resolved 
intracule densities 

^""'(u) = ^^ydR*(R)'5(r,, -u)<5,,,A,,.', (42) 

where the spin coordinates of the N — + Ni elec- 
trons in the spin- assigned wave function \E'(R) are fixed 
at Si =1 for i = 1,...,N^ and Si =| for i = + 
l,...,iV| -f- Ni. Not surprisingly, the intracule density 
at short distance is dominated by the OS contribution, 
while at long distance the SS and OS contributions be- 
come nearly identical. 

Figure 4 shows the correlation part of the radial in- 
tracule density 47rM^ — Ihf (")] , also called corre- 
lation hole, using the ZVZB2 improved estimator and 
the same Jastrow-Slater CAS(4,4) trial wave function. 
Correlation decreases the radial probability density at 



FIG. 4: Correlation part of the radial position intracule den- 
sity Attu^ [^(w) — /hf(w)] as a function of the electron-electron 
distance u for the C atom, where I{u) has been calculated 
in VMC using the ZVZB2 (with ( = 2^21 ^ 1.82 [114]) 
improved estimator and a fully optimized Jastrow-Slater 
CAS(4,4) trial wave function. The same spin contribution 
u) — /hf(w)] and the opposite spin contribution 
47rM^ [-f°^(u) - /§|(u)] are also shown. 
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FIG. 5: Correlation part of the radial position intracule den- 
sity Attu^ [^(w) — IiiF{u)] as a function of the electron-electron 
distance u for the C atom, where I{u) has been calculated 
in VMC using the ZVZB2 (with C = 2^2/ ^ 1.82 [114]) im- 
proved estimator for a series of trial wave functions of increas- 
ing accuracy: HF, Jastrow x HF, fully-optimized Jastrow x 
SD and Jastrow x CAS(4,4) wave functions. 



the two previously-mentioned maxima and increases 
it at longer distances. The same spin contribution 
Attu'^ [I^^{u) — /hf(")] ^'^d the opposite spin contribu- 
tion 47ru^ [/°^(u) — /§p(u)] are also shown in Fig. 4. 
The OS component constitutes the main contribution to 
the correlation hole. 

The effect of the chosen trial wave function on the ac- 
curacy on the correlation hole is examined in Fig. 5. Four 
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FIG. 6: Integrand of the electron-electron Coulomb interac- 
tion energy Anuliu) versus the electron-electron distance u 
for the C2 molecule, calculated in VMC using the ZVZB2 
(with C, = 2\pli ~ 1.83 [114]) improved estimator and a fully 
optimized Jastrow-Slater CAS(8,8) trial wave function. The 
same spin contribution 47rM7^^('u) and the opposite spin con- 
tribution 47rit/°^(u) are also shown. 
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FIG. 7: Correlation part of the radial position intracule den- 
sity 47ru^ [^(w) — /hf(m)] as a function of the electron-electron 
distance u for the C2 molecule, where I(u) has been cal- 
culated in VMC using the ZV2 improved estimator (with 
C = 2\/27 « 1.83 [114]) with a MCSCF CAS(8,8) wave func- 
tion and a series of Jastrow x CAS(8,8) wave functions with 
different levels of optimization. 




trial wave functions of increasing accuracy have been 
tested: HF, Jastrow x HF (with only Jastrow parameters 
optimized), Jastrow x single-determinant (SD) (with 
Jastrow and orbital parameters optimized) and Jastrow 
X CAS (4,4) (with Jastrow, CSF and orbital parameters 
optimized). Remarkably, the ZVZB2 estimator gives a 
correlation hole with a correct overall structure even with 
the uncorrelated HF wave function. For more quantita- 
tive predictions, the use of a Jastrow factor is however 
necessary. Optimization of the orbitals in a Jastrow- 
Slater single-determinant wave function brings a further 
significant improvement in the interesting short-range 
region (u < 2). The Jastrow-Slater multi-determinant 
CAS (4,4) and the Jastrow-Slater single-determinant in- 
tracules agree closely. 

C. C2 molecule 

We now discuss the more difficult case of the C2 
molecule. The ground-state wave function of this system 
has a strong multi-configurational character due to ener- 
getic near-degeneracies among the valence orbitals ( "non- 
dynamical correlation" in chemists' jargon, or "strong 
correlation" in physicists' jargon), making it a challeng- 
ing system despite its small size. 

Figure 6 plots ^-kuKu) versus u using the ZVZB2 esti- 
mator with an accurate Jastrow-Slater CAS(8,8) wave 
function (Jastrow, CSF and orbital parameters opti- 
mized). The curve displays three maxima: the maximum 
at short distance [u ~ 0.2) corresponds essentially to the 
intra-atomic K-K electron pairs; the maximum at long 
distance (u « 2.3) is located at about the interatomic 
distance and corresponds mainly to the interatomic K- 



K electron pairs; the maximum at intermediate distance 
(u ~ 1.2) must correspond mainly to intra-atomic and in- 
teratomic K-L and L-L electron pairs. As for the C atom, 
the intracule density at short distance is dominated by 
the OS contribution, while at long distance the SS and 
OS contributions become nearly identical. 

Figure 7 shows the correlation hole of the C2 molecule 
using the ZV2 improved estimator calculated in VMC 
with a MCSCF CAS(8,8) wave function and a series of 
Jastrow X CAS(8,8) wave functions for different levels of 
optimization. The MCSCF CAS (8,8) wave function does 
not correlate the core electrons and thus gives essentially 
no correlation hole at distances corresponding to intra- 
atomic core electron pairs (u ~ 0.2). At longer valence 
distances, the MCSCF CAS(8,8) wave function gives only 
a gross overall shape of the correlation hole. Introduction 
of a Jastrow factor yields a correct correlation hole at core 
distances, which as expected is about twice as deep as 
the core correlation hole of the C atom, and reduces the 
correlation hole at valence distances. The action of the 
Jastrow factor is thus not limited to very short electron- 
electron distances, but in fact importantly modifies the 
correlation hole up to distances m w 4. The Jastrow fac- 
tor has also an indirect action via reoptimization of the 
CSF and orbital coefficients in its presence. The reop- 
timization of the CSF coefficients changes slightly the 
correlation hole at valence distances. The reoptimization 
of the CSF and orbital coefficients has a more important 
impact of the correlation hole at valence distances and 
also at core distances. 

We next scrutinize in detail the convergence of the in- 
tracule density with respect to the determinantal part of 
the trial wave function. Figure 8a shows the correlation 
hole using the ZVZB2 improved estimator calculated in 
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Correlation part of the radial position intracule density Attu^ [^(w) — Ihf{u)] as a function of the electron-electron 
u for the C2 molecule, where I{u) has been calculated in VMC (a) and FN-DMC (b) using the ZVZB2 (with ( = 
1.83 [114]) improved estimator for a series of trial wave functions of increasing accuracy: Jastrow x HF, fully-optimized 
X SD, Jastrow x CAS(8,5), Jastrow x CAS(8,7) and Jastrow x CAS(8,8) wave functions. 



VMC with five trial Jastrow-Slater wave functions of in- 
creasing accuracy: Jastrow x HF (with only the Jas- 
trow parameters optimized), Jastrow x SD (with the 
Jastrow and orbital parameters optimized), and multi- 
determinant Jastrow X CAS(8,5), Jastrow x CAS(8,7) 
and Jastrow x CAS(8,8) (with the Jastrow, CSF and or- 
bital parameters optimized). At core distances {u « 0.2), 
the correlation hole is essentially converged with only 
a single-determinant Jastrow-Slater wave function, pro- 
vided that the orbitals are reoptimized together with the 
Jastrow factor. In contrast, at longer valence distances, 
the correlation hole depends very strongly on the deter- 
minantal part of the wave function. In particular, the 
multi-determinant wave functions yield a depletion of the 
intracule density at about the bond length, a feature not 
present at the single-determinant level. We have verified 
that this minimum disappears when using a pseudopo- 
tential removing the Is electrons, and we thus interpret 
it as a decrease of probability of finding two interatomic 
core electrons separated by the bond distance. Overall, 
it appears necessary to use at least a multi-determinant 
CAS (8, 7) wave function which includes configurations 
constructed from the antibonding tt orbitals to reach rea- 
sonable convergence. We have checked that further exci- 
tations beyond the CAS(8,8) wave function barely change 
the correlation hole. The corresponding correlation holes 
calculated in FN-DMC with the same ZVZB2 estimator 
and trial wave functions are reported in Fig. 8b. The 
DMC intracules do not differ much from the correspond- 
ing VMC intracules, the accuracy of the intracule densi- 
ties in the valence region being still essentially controlled 
by the trial wave function. 
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FIG. 9: Correlation part of the radial position intracule den- 
sity Attu^ [^(^') ~ Ihf{u)] as a function of the electron-electron 
distance u for the N2 molecule, where I{u) has been calcu- 
lated in VMC using the ZVZB2 (with ( = 2V2I a; 2.14 [114]) 
improved estimator for a series of trial wave functions of in- 
creasing accuracy: Jastrow x HF, fully-optimized Jastrow x 
SD and Jastrow x CAS(10,8) wave functions. 



D. N2 molecule 

We finish our illustration of the calculation of intrac- 
ule densities with the N2 molecule, which has recently 
been the object of experimental and theoretical investi- 
gations [42]. 

The correlation holes of this molecule calculated in 
VMC with the ZVZB2 improved estimator for a Jastrow 
x HF, and fully optimized Jastrow x SD and Jastrow x 
CAS(10,8) wave functions are plotted in Fig. 9. In sharp 
contrast to the C2 molecule, no depletion of intracule 
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density is observed at the bond length. Also, the correla- 
tion hole is here essentially converged with only a single- 
determinant Jastrow-Slater wave function provided that 
the orbitals are optimized with the Jastrow factor. Unlike 
the C2 molecule, going to a multi-determinant wave func- 
tion does not have a large effect on the correlation hole 
at valence distances. The correlation hole calculated here 
agrees reasonably well with recent multi-reference config- 
uration interaction and coupled-cluster calculations, and 
with experiment [42]. 

VI. CONCLUSIONS 

We have presented improved QMC estimators for the 
spherically averaged position intracule density I{u), con- 
structed using the general zero- variance zero-bias princi- 
ple for observables that do not commute with the Hamil- 
tonian. By replacing the average of the local delta- 
function operator by the average of a smooth non-local 
operator, these estimators decrease the variance of the 
standard histogram estimator by several orders of mag- 
nitude, and thus make the calculation of this quantity 
in QMC vastly more efficient. Interestingly, they per- 
mit calculations of I{u) for very short and very large 
intcrclectronic distances u that are never realized in the 
Monte Carlo run. These new estimators can also decrease 
the systematic error of the intracule density due to the 
approximate trial wave function. Other advantages of 
these estimators are the absence of any discretization er- 
ror with respect to u and the possibility to obtain very 
smooth curves for I{u). These improved estimators, to- 
gether with the achievement of systematically reducing 



the systematic error in both VMC and DMC calculations 
by optimization of trial wave functions with an increasing 
number of parameters, have allowed us to obtain accurate 
correlated intracule densities for atoms and molecules. 

The estimators presented here can be used with triv- 
ial adaptations for QMC calculations of the companion 
entity of I{u), namely the spherically averaged extracule 
density E{r), representing the probability density of find- 
ing two electrons with center of mass at a radial distance 
r with respect to the chosen origin. Similar improved 
estimators can be constructed for the three-dimensional 
intracule density /(u), extracule density E(r) and the 
full pair density n2(ri,r2). More generally, the variance 
reduction technique presented here can be applied to cal- 
culations of any pair-correlation function in classical and 
quantum Monte Carlo calculations. 
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